/**
 * @author	Jiefei Ma
 * @date		March 2010
 * Department of Computing, Imperial College London
 * 
 * NOTE: this version can be run on SICStus4 or YAP6
 * NOTE: this version contains two meta-interpreters -- depth-first and 
 *       heuristics search
 * NOTE: this version attempts to compete with existing abductive system
 *       on performance.
 */


/* 
 *      Added for the learning:
 *
 *      Each state carries a Info field. For the learning the field contains a list of pairs, e.g.
 *      Info = [(depth, 51), (score, 3)]
 *
 *      In case the debug mode is active info contains an additional field e.g. 
 *      Info = [(depth, 51), (score, 3), (dbg, 309-[goal])] 
 */


% # ===================== End-User Documentation ======================
% 
% 1. Syntax for user input files:
% 
% * abducible, defined (a.k.a. non-abducible) and builtin are predicates
%   e.g. holds(has(jeff, apple), 3), member(X, [2,3])
%
% * a negative literal is "\+ A" where "A" is an abducible, defined or builtin
%
% * a finite domain constraint is one used in "library(clpfd)"
%   e.g. X in 1..3, Y #< Z - 3, X #\= Y
%
% * a real domain constraint is one used in "library(clpr)"
%   e.g. {X + 3 = Y -2}
%
% * an equality is "X = Y", an inequality is "X =/= Y".
%
% 2. Each input file should specify an abductive framework:
%
% * an abducible is declared as:
%
%   abducible(happens(_,_)).
%
% * a builtin is declared as:
%
%   builtin(member(_,_)).
%
% * a rule is specified as:
%
%   holds(F,T) :- 
%     happens(A, T1),
%     T1 #< T,
%     initiates(A, F, T1),
%     \+ clipped(F, T1, T).
%
% * an integrity is specified as:
%
%   ic :- happens(A1, T), happens(A2, T), A1=/=A2.
%
% 3. Loading and querying the System
%
% * Load a user defined abductive theory of <P, A, B, I>, where "P" is a 
%   set of rules called the background knowledge, "A" is the set of abducible 
%   predicates, "B" is the set of builtin predicates, and "I" is the set of 
%   integrity constraints,  e.g. 
%     ?- load_theory('blocksworld.pl').
%
% * Query the system:
%   e.g. 
%     ?- query([T #> 0, T #< 10, holds(has(jiefei, apple), T)], Ans)
%
%  If succeeded, Ans has the form (As, (Es, Fs, Rs), Ns, Info) where "As" is 
%  the set of collected abducibles, "Es" is the set of collected inequalities,
%  "Fs" is the set of finite domain constraints, "Rs" is the set of real domain
%  constraints, "Ns" is the set of collected dynamic denials and "Info" is the 
%  user collected data along the computation.
%  
% ==========================  End of Documentation ===========================

:- multifile builtin/1.
:- discontiguous wrap_atom/2.

:- use_module(library(clpfd)).
:- use_module(library(clpr)).
:- use_module(inequalities).
% For solver hooks
:- use_module(library(atts)).
% For optimisations
:- use_module(library(ordsets), [
                list_to_ord_set/2
        ]).

:- [heuristics].

% ===================== End-User API ================================
load_theory(DataFile) :-
	loadf(DataFile, Cls),
	assert_knowledge(Cls).

load_theories(DataFiles) :-
	loadfs(DataFiles, Cls),
	assert_knowledge(Cls).

clear_theory :-
	clean_up.

% Ans: (As, Cs, Ns, Info)
query(Query, Ans) :-
	initialise(Query, Gs, As, Cs, Ns, Info),
	solve_all(Gs, As, Cs, Ns, Info, Ans).

% Ans: Score-Query-(As, Cs, Ns, Info)
heuristics_query(Query, Ans) :-
	initialise(Query, Gs, As, Cs, Ns, Info),
	heuristics_reset,
        init_learning_score(InitScore),
        set_info(Info, Info2, scr, InitScore),
        InitState = (Gs, As, Cs, Ns, Query, Info2),
	new_state_candidate_store(EmptyStore),
	add_states([InitScore-InitState], EmptyStore, InitStore),
	h_solve_all(InitStore, Ans).

debug_on :-
	set_value(dbg, true),
	set_value(max_state, 0).

debug_off :-
	set_value(dbg, false).

debug_set_max_state(X) :-
	integer(X),
	X >= 0,
	set_value(max_state, X).

debug_trace(File) :-
	open(File, write, OUT),
	catch((construct_and_dump_tree(OUT)), Err, (write(Err), nl)),
	flush_output(OUT),
	close(OUT).

enforce_labeling(true) :-
	set_value(lbl, true).
enforce_labeling(false) :-
	set_value(lbl, false).

partial_evaluation(true) :-
	set_value(prep, true).
partial_evaluation(false) :-
	set_value(prep, false).

partial_evaluation_set_max_depth(X) :-
	integer(X),
	X > 0,
	set_value(prep_max_depth, X).

% ===================== Implementation ==============================

% ### Preamable ###





% need to declare the following operator for the subsequent conditional compilations
:- if(current_prolog_flag(dialect, yap)).
:- op(760, yfx, #<=>).
:- elif(current_prolog_flag(dialect, sicstus)).
:- op(760, yfx, #<==>).
:- endif.

:- if(current_prolog_flag(dialect, yap)).
%{
fd_entailment(C, B) :- 
	call(C #<==> B).

finite_domain_inequality(X=/=Y) :-
	fd_var(X) ; fd_var(Y) ; % already be fd var
	(var(X), integer(Y)) ; (var(Y), integer(X)). % about to become fd var

real_domain_inequality(X=/=Y) :-
	get_attr(X, itf, _) ; get_attr(Y, itf, _) ; % already be rd var
	(var(X), float(Y)) ; (var(Y), float(X)). % about to become rd var

%}
:- elif(current_prolog_flag(dialect, sicstus)).
%{
fd_entailment(C, B) :-
	call(C #<=> B).

finite_domain_inequality(X=/=Y) :-
	fd_var(X) ; fd_var(Y) ; % already be fd var
	(var(X), integer(Y)) ; (var(Y), integer(X)). % about to become fd var

real_domain_inequality(X=/=Y) :-
	% get_atts(X, itf) ; get_atts(Y, itf) ; % already be rd var % FIXME: can't get it working
	(var(X), float(Y)) ; (var(Y), float(X)). % about to become rd var

%}
:- endif.

% *** Utilities ***
:- if(current_prolog_flag(dialect, yap)).
%{
:- use_module(library(terms), [
		unifiable/3,
		new_variables_in_term/3,
		variables_within_term/3
	]).
:- use_module(library(apply_macros), [
		maplist/3,
		selectlist/3
	]).
:- use_module(library(lists), [
		member/2,
		append/3,
		select/3,
		reverse/2
	]).
rev(X, Y) :- reverse(X, Y).
%}
:- elif(current_prolog_flag(dialect, sicstus)).
%{
:- use_module(library(terms), [
		term_variables/2
	]).
:- use_module(library(lists), [
		maplist/3,
		select/3,
		rev/2
	]).

variables_within_term(Vs, Term, OldVs) :-
	term_variables(Term, TVs),
	collect_old_variables(TVs, Vs, OldVs).
collect_old_variables([], _, []).
collect_old_variables([H|T], Vs, OldVs) :-
	collect_old_variables(T, Vs, OldVs1),
	(strictmember(Vs, H) ->
		OldVs = [H|OldVs1]
	;
		OldVs = OldVs1
	).

unifiable(X, Y, Eq) :-
	(var(X) ; var(Y)), !,
	(X == Y -> Eq = [] ; Eq = [X = Y]). 
unifiable(X, Y, []) :-
	atomic(X), atomic(Y), !, X == Y.
unifiable(X, Y, Eqs) :-
	functor(X, F, A),
	functor(Y, F, A), 
	X =.. [F|ArgsX],
	Y =.. [F|ArgsY],
	all_unifiable(ArgsX, ArgsY, Eqs).
all_unifiable([], [], []).
all_unifiable([X|TX], [Y|TY], AllEqs) :-
	unifiable(X, Y, Eqs),
	all_unifiable(TX, TY, RestEqs),
	append(Eqs, RestEqs, AllEqs).

selectlist(Pred, L1, L2) :- 
	selectlist_aux(L1, Pred, L2).
selectlist_aux([], _, []).
selectlist_aux([H|T], P, L) :-
	selectlist_aux(T, P, L1),
	(\+ \+ call(P, H) ->
		L = [H|L1]
	;
		L = L1
	).

:- dynamic mysetvalue/2.
set_value(Key, Val) :- retractall(mysetvalue(Key, _Val)), assert(mysetvalue(Key, Val)).
get_value(Key, Val) :- mysetvalue(Key, Val).

forall(X, Y) :- \+ (X, \+ Y).
%}
:- endif.

unifiable(X, Y) :- \+ \+ (X = Y).

nonground(X) :- \+ ground(X).

strictdelete([], _, []).
strictdelete([H|T], X, R) :-
	(H == X ->
		strictdelete(T, X, R)
	;
		strictdelete(T, X, R1),
		R = [H|R1]
	).

strictmember([H|T], X) :-
	(X == H ->
		true
	;
		strictmember(T, X)
	).

% ----------------------------------------------------------------------------------------
% --- Preprocessing ---

% Predicates used for specifying the in
:- dynamic abducible/1, builtin/1, control/1. % used by users and the system
:- dynamic rule/2, ic/1. % used by the meta-interpreter
:- dynamic solution_state/1.

clean_up :-
	retractall(abducible(_)),
	retractall(builtin(_)),
        retractall(control(_)),
	retractall(rule(_,_)),
	retractall(ic(_)).

% [Literal Type]: 
%
% In order to speed up the search, each input rule or integrity constraint
% is preprocessed when it is loaded, such that each literal "L" is wrapped 
% as "(Type, L)", where "Type" is the literal type.
%
% There are four literal types: 
%  a: abducible
%  b: builtin
%  d: defined
%  e: equality
%  i: inequality
%  n: negation
%  f: finite domain constraint
%  r: real domain constraint
%  -: failure goal (see below)
%
% There is a special literal, of the form "fail(Vs, Literals)", called the 
% failure goal which can appear as a sub-goal during the abductive search.  
% Logically, the failure goal is "\forall Vs . <- Literals", where "Vs" is 
% the set of variables in "Literals" that are universally quantified with
% the scope the whole failure goal (denial).  All other varaibles in "Literals"
% are existentially quantified implicitly. "-" represents the type of a failure
% goal.

wrap_atom(X, (b, X)) :-
        builtin(X), !.

wrap_literal(X, Xw) :-
	(X = (\+ A) ->
		% for negative literal, we need to wrap the atom too
		wrap_atom(A, Aw),
		Xw = (n, \+ Aw)
	;
		wrap_atom(X, Xw)
	).

wrap_atom(X = Y, (e, X = Y)) :- !.
wrap_atom(X=/=Y, (i, X=/=Y)) :- !.

wrap_atom({X}, (r, {X})) :- !.

wrap_atom(X in Y, (f, X in Y)) :- !.
wrap_atom(X #= Y, (f, X #= Y)) :- !.
wrap_atom(X #< Y, (f, X #< Y)) :- !.
wrap_atom(X #=< Y, (f, X #=< Y)) :- !.
wrap_atom(X #> Y, (f, X #> Y)) :- !.
wrap_atom(X #>= Y, (f, X #>= Y)) :- !.
wrap_atom(X #\= Y, (f, X #\= Y)) :- !.
wrap_atom(#\ X, (f, #\ X)) :- !.
wrap_atom(X #/\ Y, (f, X #/\ Y)) :- !.
wrap_atom(X #\/ Y, (f, X #\/ Y)) :- !.

wrap_atom(X, (a, X)) :- 
	abducible(X), !.
wrap_atom(X, (d, X)). % for anything else, it is assumed to be a defined.

unwrap_literal((n, \+ (_, A)), (\+ A)) :- !.
unwrap_literal((-, fail(Vs, BodyW)), fail(Vs, Body)) :- 
	!, maplist(unwrap_literal, BodyW, Body).
unwrap_literal((_, A), A).

flatten_body(B, L) :-
	(B = (X, Y) ->
		flatten_body(Y, L1),
		L = [X|L1]
	;
		L = [B]
	).

wrap_ic(Body, ic(NewBody)) :-
	flatten_body(Body, Lits),
	maplist(wrap_literal, Lits, NewBody).
wrap_fact(Head, rule(NewHead, [])) :-
	wrap_atom(Head, NewHead).
wrap_rule(Head, Body, rule(NewHead, NewBody)) :-
	flatten_body(Body, Lits),
	maplist(wrap_literal, [Head|Lits], [NewHead|NewBody]).

unwrap_denial(fail(Vs, BodyW), fail(Vars, Body)) :-
	maplist(unwrap_literal, BodyW, Body),
	selectlist(nonground, Vs, Vars).
		
% /** API **/
% Read the user input file and perform preprocessing.
loadfs([], []).
loadfs([F|T], Cls):-
	loadfs(T, C2),
	loadf(F, C1),
	append(C1, C2, Cls).

loadf(DataFile, Cls) :-
	open(DataFile, read, IN),
	read_clauses(IN, Cls),
	close(IN).

assert_knowledge(Cls) :-
	assert_abducibles(Cls, Cls1),
	assert_builtins(Cls1, Cls2),
        assert_controls(Cls2, Cls3),
	assert_clauses(Cls3),
	(get_value(prep, true) ->
		% allows partial evaluation of the theory
		findall(rule(H, B), rule(H, B), Rules),
		partial_evaluate(Rules, [], NewRules, _),
		% replace the old theory with the equivalent but new theory
		retractall(rule(_, _)),
		forall(member(R, NewRules), assert(R))
	;
		true
	).

read_clauses(IN, Cls) :-
	catch(read(IN, Term), Err, (write(Err), fail)),
	(Term == end_of_file ->
		Cls = []
	; Term = (:- _D) ->
		%call(D),
		read_clauses(IN, Cls)
	;
		read_clauses(IN, ClsRest),
		Cls = [Term|ClsRest]
	).

assert_abducibles([], []).
assert_abducibles([abducible(X)|T], L) :-
	!, assertz(abducible(X)),
	assert_abducibles(T, L).
assert_abducibles([H|T], [H|L]) :-
	assert_abducibles(T, L).

assert_builtins([], []).
assert_builtins([builtin(X)|T], L) :-
	!, assertz(builtin(X)),
	assert_builtins(T, L).
assert_builtins([H|T], [H|L]) :-
	assert_builtins(T, L).

assert_controls([], []).
assert_controls([control(X)|T], L) :-
        !, assertz(control(X)),
        assert_controls(T, L).
assert_controls([H|T], [H|L]) :-
        assert_controls(T, L).

assert_clauses([]).
assert_clauses([C|T]) :-
	(C = ( H :- B) ->
		(H == ic ->
			% preprocess and add integrity constraint
			wrap_ic(B, IC),
			assert(IC)
		;
			% preprocess and add a rule
			wrap_rule(H, B, Rule),
			assert(Rule)
		)
	;
		wrap_fact(C, Fact),
		assert(Fact)
	),
	assert_clauses(T).

% =====================  Debug Facilities ===========================
:- dynamic arc/2, state/1, solved/1, failed/1.

reset_state_counter :-
	set_value(counter, 0).
new_state_id(ID) :-
	get_value(counter, ID),
	NextID is ID + 1, 
	set_value(counter, NextID).

debug_hook(init_state(Query), Info) :-
	!, 
	(get_value(dbg, true) ->
		% clean up and reset
		retractall(arc(_,_)),
		retractall(state(_)),
		retractall(solved(_)),
		reset_state_counter,
		% record infor
		new_state_id(RootID),
		RootID = 0,
		set_info([], Info, dbg, RootID-Query)
	;
		Info = []
	).
debug_hook(inter_state(OldInfo, State), NewInfo) :-
	!,
	(get_value(dbg, true) ->
		get_info(OldInfo, dbg, ParentID-Query),
		new_state_id(ID),
                %(ParentID = 677 -> trace ; true),
		((get_value(max_state, Max), Max > 0, ID > Max) ->
			abort
		;
			true
		),
		assert(arc(ParentID, ID)),
		State = (Gs, As, Cs, Ns),
		extract_constraints(Cs, NewCs),
		NewState = (Gs, As, NewCs, Ns),
		assert(state(ID-Query-NewState-OldInfo)),
		set_info(OldInfo, NewInfo, dbg, ID-Query)
	;
		NewInfo = OldInfo
	).
debug_hook(solved_state(OldInfo, State), NewInfo) :-
	(get_value(dbg, true) ->
		get_info(OldInfo, dbg, ParentID-Query),
		new_state_id(ID),
%(ParentID = 1258 -> trace ; true),
                assert(arc(ParentID, ID)),
		assert(solved(ID)),
		State = (Gs, As, Cs, Ns),
		extract_constraints(Cs, NewCs),
		NewState = (Gs, As, NewCs, Ns),
		assert(state(ID-Query-NewState-OldInfo)),
		set_info(OldInfo, NewInfo, dbg, ID-Query)
	;
		NewInfo = OldInfo
	).
construct_and_dump_tree(OUT) :-
	write(OUT, ':- use_module(library(clpfd)).'), nl(OUT),
	write(OUT, ':- use_module(library(clpr)).'), nl(OUT),
	write(OUT, ':- use_module(inequalities).'), nl(OUT),
	write(OUT, '% arc(ParentStateID, ChildStateID).'), nl(OUT), 
	forall(arc(X,Y), (X \= 0 -> portray_clause(OUT, arc(X,Y)); true)),
	write(OUT, '% st(StateID, StateType, Query-(Goals, Abducibles, Constraints, Denials), Info)'), nl(OUT),
	write(OUT, '%   StateType: i (intermediate), f (failed), s (succeeded)'), nl(OUT),
	write(OUT, '%   Info: a list of key-value pair.  In this version, it is always empty.'), nl(OUT),
	forall(state(ID-Q-State-Info), (
		State = (Gw, A, Cs, Nw),
		maplist(unwrap_literal, Gw, G),
		maplist(unwrap_denial, Nw, N),
		(solved(ID) ->
			Type = s
		; \+ arc(ID, _)->
			Type = f
		;
			Type = i
		),
		portray_clause(OUT, st(ID, Type, Q-(G, A, Cs, N), Info))
	)).

% --- Abductive Meta-Interpreter (Depth-First) ---

% solve_all(Gs, As, Cs, Ns, Info, Solution)
%   Gs: set of pending goals
%   As: set of collected abducibles
%   Cs: set of constraints, of the form (Es, Fs, Rs), where
%     Es: set of inequalities
%     Fs: set of finite domain constraints
%     Rs: set of real domain constraints
%   Info: any user-defined data structure for storing 
%     computational information, such as debug data
%   Solution: when the whole abductive computation succeeds,
%     this will contain the final answer

% base case: no more pending goals.  so it succeeds.
solve_all([], As, (Es, Fs, Rs), NsW, Info, (As, Cs, Ns, NewInfo)) :-
	debug_hook(solved_state(Info, ([], As, (Es, Fs, Rs), NsW)), NewInfo),
	extract_constraints((Es, Fs, Rs), Cs),
	maplist(unwrap_denial, NsW, Ns).
% recursive case: simple depth-first left-right search strategy
solve_all([(Type, G)|RestGs], As, Cs, Ns, Info, Sol) :-
	((Type \= -) ->
		debug_hook(inter_state(Info, ([(Type, G)|RestGs], As, Cs, Ns)), NewInfo) 
	; 
		NewInfo = Info
	),
	solve_one(Type, G, RestGs, As, Cs, Ns, NewInfo, Sol).

solve_one(a, A, RestGs, As, Cs, Ns, Info, Sol) :-
	(
		% reuse hypotheses
		member(A, As),
		solve_all(RestGs, As, Cs, Ns, Info, Sol)
	;	% a create new hypothesis
		resolve_abducible_with_delta(As, A, Inequalities), 
		propagate_inequalities(Inequalities, Cs, NewCs), 
		resolve_abducible_with_dynamic_ics(Ns, A, FailureGoals), 
		append(FailureGoals, RestGs, NewGs),
		solve_all(NewGs, [A|As], NewCs, Ns, Info, Sol)
	).

solve_one(d, D, RestGs, As, Cs, Ns, Info, Sol) :-
	rule((d,D), B), % pick a rule
	append(B, RestGs, NewGs), % FIXME: solve constraints first?
	solve_all(NewGs, As, Cs, Ns, Info, Sol).

solve_one(b, B, RestGs, As, Cs, Ns, Info, Sol) :-
	call(B), % backtrackable
	solve_all(RestGs, As, Cs, Ns, Info, Sol).

solve_one(e, X = Y, RestGs, As, Cs, Ns, Info, Sol) :-
	call(X = Y),
	solve_all(RestGs, As, Cs, Ns, Info, Sol).

solve_one(i, X=/=Y, RestGs, As, Cs, Ns, Info, Sol) :-
	propagate_inequalities([X=/=Y], Cs, NewCs), 
	solve_all(RestGs, As, NewCs, Ns, Info, Sol).

solve_one(n, \+ G, RestGs, As, Cs, Ns, Info, Sol) :-
	solve_one(-, fail([], [G]), RestGs, As, Cs, Ns, Info, Sol).

solve_one(f, FC, RestGs, As, Cs, Ns, Info, Sol) :-
	propagate_finite_domain_constraints([FC], Cs, NewCs),
	solve_all(RestGs, As, NewCs, Ns, Info, Sol).

solve_one(r, RC, RestGs, As, Cs, Ns, Info, Sol) :-
	propagate_real_domain_constraints([RC], Cs, NewCs), 
	solve_all(RestGs, As, NewCs, Ns, Info, Sol).

solve_one(-, fail(Vs, Lits), RestGs, As, Cs, Ns, Info, Sol) :-
	debug_hook(inter_state(Info, ([(-, fail(Vs, Lits))|RestGs], As, Cs, Ns)), NewInfo), 
	(safe_select_failure_literal(Lits, Vs, (Type, L), RestLits) ->
		fail_one(Type, L, Vs, RestLits, RestGs, As, Cs, Ns, NewInfo, Sol)
	;
		(Lits == [] ->
			fail % fail normally
		;
			% if "Lits" contains only finite domain constraints and/or 
			% real constraints, then we can try to test for satisfiability.
			% if it is satisfiable, then it fails normally (i.e. falsity 
			% can be derived); otherwise if it is not satisfiable, then
			% it succeeds.
			(partition_failure_literals(Lits, FCs, RCs) ->
				(\+ \+ (
							% can be satisfied?
							propagate_finite_domain_constraints(FCs, Cs, Cs1),
							propagate_real_domain_constraints(RCs, Cs1, _)
						) ->
					% yes, satisfiable.  so fail
					fail
				;
					% no, and hence no floundering, so continue the reasoning
					solve_all(RestGs, As, Cs, Ns, NewInfo, Sol)
				)
			; % indeed, it flounders.
				write('Floundered: '), write(fail(Vs, Lits)), nl,
				fail
			)
		)
	).

fail_one(a, A, Vs, RestLits, RestGs, As, Cs, Ns, Info, Sol) :-
	resolve_failure_abducible_with_delta(As, A, Vs, RestLits, FailureGoals),
	append(FailureGoals, RestGs, NewGs),
	NewNs = [fail(Vs, [(a, A)|RestLits])|Ns],
	solve_all(NewGs, As, Cs, NewNs, Info, Sol).

fail_one(d, D, Vs, RestLits, RestGs, As, Cs, Ns, Info, Sol) :-
	findall(H-B, (rule((d, H), B), unifiable(H, D)), Rules),
	resolve_failure_non_abducible_with_rules(Rules, D, Vs, RestLits, FailureGoals),
	append(FailureGoals, RestGs, NewGs),
	solve_all(NewGs, As, Cs, Ns, Info, Sol).

fail_one(b, B, Vs, RestLits, RestGs, As, Cs, Ns, Info, Sol) :-
	findall(B-[], call(B), Rules),
	resolve_failure_non_abducible_with_rules(Rules, B, Vs, RestLits, FailureGoals),
	append(FailureGoals, RestGs, NewGs),
	solve_all(NewGs, As, Cs, Ns, Info, Sol).

fail_one(e, X = Y, Vs, RestLits, RestGs, As, Cs, Ns, Info, Sol) :-
	((var(X), strictmember(Vs, X)) ->
		% X is uni. quant. and we don't care about Y
		strictdelete(Vs, X, NewVs),
		call(X = Y),
		solve_one(-, fail(NewVs, RestLits), RestGs, As, Cs, Ns, Info, Sol)
	; (var(Y), strictmember(Vs, Y)) ->
		% Y is uni. quant. but X is not
		strictdelete(Vs, Y, NewVs),
		call(Y = X),
		solve_one(-, fail(NewVs, RestLits), RestGs, As, Cs, Ns, Info, Sol)
	; var(X) ->
		% X is ex. quant. 
		% NB: by the safe selection strategy, "Y" doesn't contain any
		% universally quantified variable
		(
			% try to succeed in the inequality
			propagate_inequalities([X=/=Y], Cs, NewCs),
			solve_all(RestGs, As, NewCs, Ns, Info, Sol)
		;
			% try to succeed in the equality and fail one of the rest literals
			call(X = Y),
			solve_one(-, fail(Vs, RestLits), RestGs, As, Cs, Ns, Info, Sol)
		)
	; var(Y) ->
		% Y is ex. quant. but X is not
		% NB: by the safe selection strategy, "X" doesn't contain any
		% universally quantified variable
		(
			% try to succeed in the inequality
			propagate_inequalities([Y=/=X], Cs, NewCs),
			solve_all(RestGs, As, NewCs, Ns, Info, Sol)
		;
			% try to succeed in the equality and fail one of the rest literals
			call(Y = X),
			solve_one(-, fail(Vs, RestLits), RestGs, As, Cs, Ns, Info, Sol)
		)
	; % one of them is a variable
		(unifiable(X, Y, Es) ->
			maplist(wrap_atom, Es, EsW),
			append(EsW, RestLits, NewBody),
			solve_one(-, fail(Vs, NewBody), RestGs, As, Cs, Ns, Info, Sol)
		;
			% this literal fails trivially
			solve_all(RestGs, As, Cs, Ns, Info, Sol)
		)
	).

fail_one(i, X=/=Y, Vs, RestLits, RestGs, As, Cs, Ns, Info, Sol) :-
	(
		call(X = Y),
		solve_all(RestGs, As, Cs, Ns, Info, Sol)
	;
		propagate_inequalities([X=/=Y], Cs, NewCs),
		solve_one(-, fail(Vs, RestLits), RestGs, As, NewCs, Ns, Info, Sol)
	).

fail_one(n, \+ G, Vs, RestLits, RestGs, As, Cs, Ns, Info, Sol) :-
	(
		% Case 1:
		solve_all([G|RestGs], As, Cs, Ns, Info, Sol)
	;
		% case 2:
                RestLits \= [], %DOMENICO-5-Sept
		G = (Type, L),
		fail_one(Type, L, [], [], [(-, fail(Vs, RestLits))|RestGs], As, Cs, Ns, Info, Sol)
	).

fail_one(f, C, Vs, RestLits, RestGs, As, Cs, Ns, Info, Sol) :-
	fd_entailment(C, B), % reitification
	(B == 0 -> 
		% the constraint can never hold anyway
		solve_all(RestGs, As, Cs, Ns, Info, Sol)
	; B == 1 ->
		% the constraint always hold
		solve_one(-, fail(Vs, RestLits), RestGs, As, Cs, Ns, Info, Sol)
	; % the constraint can hold and can unhold
		% then we either (a) fail it or (b) succeed it and fail the rest
		(
			negate_finite_domain_constraint(C, NC), 
			propagate_finite_domain_constraints([NC], Cs, NewCs),
			solve_all(RestGs, As, NewCs, Ns, Info, Sol)
		;
			propagate_finite_domain_constraints([C], Cs, NewCs),
			solve_one(-, fail(Vs, RestLits), RestGs, As, NewCs, Ns, Info, Sol)
		)
	).

fail_one(r, C, Vs, RestLits, RestGs, As, Cs, Ns, Info, Sol) :-
	negate_real_domain_constraint(C, NC),
	C = {Cons}, NC = {NCons},
	(entailed(NCons) -> 
		% the constraint can never hold anyway
		solve_all(RestGs, As, Cs, Ns, Info, Sol)
	; entailed(Cons) ->
		% the constraint always hold
		solve_one(-, fail(Vs, RestLits), RestGs, As, Cs, Ns, Info, Sol)
	; % we either (a) fail it, or (b) succeed it and fail the rest.
		(
			propagate_real_domain_constraints([NC], Cs, NewCs),
			solve_all(RestGs, As, NewCs, Ns, Info, Sol)
		;
			propagate_real_domain_constraints([C], Cs, NewCs),
			solve_one(-, fail(Vs, RestLits), RestGs, As, NewCs, Ns, Info, Sol)
		)
	).

% --- Abductive Meta-Interpreter (Heuristic Search) ---


display_stats_in_derivation(ASL, ANSN, L) :-
        get_stat(explored_states, L1),
        print('\n###########################'),
        print('\n# Current situation'),
        print('\n# Open states: '), print(ASL),
        print('\n# Open states on average: '), print(ANSN),
        print('\n# Iterations : '), print(L),
        print('\n# Total explored states : '), print(L1),
        print('\n###########################\n').


% State = (Gs, As, Cs, Ns, Q, Info)
%   Gs: set of pending goals
%   As: set of collected abducibles
%   Cs: set of constraints, of the form (Es, Fs, Rs), where
%     Es: set of inequalities
%     Fs: set of finite domain constraints
%     Rs: set of real domain constraints
%   Q: original query
%   Info: any user-defined data structure for storing 
%     computational information, such as debug data

h_solve_all(_States, Solution) :-
        terminate(Solution), !.
        
h_solve_all(States, Solution) :-
	select_state(States, Score-State, RestStates),
	State = (Gs, As, Cs, Ns, Q, Info),
%        portray_list(As), nl, nl,
	(inject_constraints(Cs) -> % HACK: to cope with unsound CLP(FD)
		% HACK: to cope with findall
		(Gs = [(Type, G)|RestGs] ->
			%print((Type, G)),
                        debug_hook(inter_state(Info, (Gs, As, Cs, Ns)), NewInfo),    %debug info updated
                        update_info_post_debug(NewInfo, NewInfo2),             %other info updated. Currently it updates the depth (for the children)
			findall((NextScore, Depth)-NextStateWithScoreUpdated, (
					h_solve_one(Type, G, RestGs, As, Cs, Ns, Q, NewInfo2, NewState),
                                        get_info(NewInfo2, depth, Depth), 
                                        %PRUNES if the depth exceeds the limits
                                        check_depth(Depth),                                                     
					update_inequalities(NewState, NextState), % HACK: to cope with findall
					copy_term(State, CState), copy_term(NextState, CNextState),
                                        copy_term(G, CG),
                                        prune((Type, CG), CState, CNextState, NextScore),
                                        learning_score((Type, CG), CState, CNextState, NextScore),
                                        update_new_state_with_new_score(NextScore, NextState, NextStateWithScoreUpdated, _UpdatedInfo)                                         
				), NewStates), 
                        add_states(NewStates, RestStates, AllStates),
                        %STATS%%%%%%%%%
                                                        length(AllStates, ASL),                                                        
                                                        length(NewStates, L),
                                                        increment_stat(explored_states, L),
                                                        (L = 0 -> increment_stat(failed_branches) ; true),
                                                        get_stat(iterations, ITS),
                                                        increment_stat(iterations),
                                                        get_stat(average_n_states, ANS),
                                                        ANSN is (ANS * ITS + ASL)/(ITS + 1),
                                                        set_stat(average_n_states,ANSN),
                                                        DisplayStat is ITS mod 500,
                                                        ((DisplayStat == 0) -> display_stats_in_derivation(ASL, ANSN, ITS) ; true),
			%%%%%%%%%%%%%%%%%
%                        findall(  (NextScore1, Depth1)-Infoo2, member( (NextScore1, Depth1)-(_Gs, _As, _Cs, _Ns, _Q, Infoo2), AllStates), QQQQ),
%                        portray_list(QQQQ), nl, nl, nl,
			h_solve_all(AllStates, Solution)
		;
			debug_hook(solved_state(Info, (Gs, As, Cs, Ns)), NewInfo),
                        update_info_post_debug(NewInfo, NewInfo2),
			( % return this
                                print('\nSolution leaf reached.\n'),
				extract_constraints(Cs, NewCs),
				maplist(unwrap_denial, Ns, NewNs),
				(Solution = Score-Q-(As, NewCs, NewNs, NewInfo2), ! ; true),                                
                                assert(solution_state(Solution)),
			%; % or find the next one
				h_solve_all(RestStates, Solution)
			)
		)
	;
		h_solve_all(RestStates, Solution)
	).

% HACK: start
inject_constraints((Es, Fs, Rs)) :-
	inject_all(Es),
	inject_all(Fs),
	inject_all(Rs).
inject_all([]).
inject_all([H|T]) :-
	call(H),
	inject_all(T).

update_inequalities(
		(Gs, As, (Es, Fs, Rs), Ns, Q, Info), 
		(Gs, As, (NewEs, Fs, Rs), Ns, Q, Info)) :-
	term_variables(Es, EsVs),
	maplist(inequalities, EsVs, AllEs),
	flatten_lists(AllEs, EsTmp),
	list_to_ord_set(EsTmp, NewEs).
% HACK: end


h_solve_one(a, arwr(A), RestGs, As, Cs, Ns, Q, Info, NextState) :- !,
        ( 
                resolve_abducible_with_dynamic_ics(Ns, A, FailureGoals), 
                sort_failure_goals(FailureGoals, SortedFailureGoals),
                append(SortedFailureGoals, RestGs, NewGs),
                NextState = (NewGs, As, Cs, Ns, Q, Info)
        ).


h_solve_one(a, A, RestGs, As, Cs, Ns, Q, Info, NextState) :-
	(
		% reuse hypotheses
		member(A, As),
		NextState = (RestGs, As, Cs, Ns, Q, Info)
	;	% a create new hypothesis
		resolve_abducible_with_delta(As, A, Inequalities), 
		propagate_inequalities(Inequalities, Cs, NewCs), 
		resolve_abducible_with_dynamic_ics(Ns, A, FailureGoals), 
                sort_failure_goals(FailureGoals, SortedFailureGoals),
		append(SortedFailureGoals, RestGs, NewGs),
		NextState = (NewGs, [A|As], NewCs, Ns, Q, Info),
                (ground([A|As]) ->
                 checkAs([A|As])
                        ;
                true)
	).



h_solve_one(d, D, RestGs, As, Cs, Ns, Q, Info, NextState) :-
	rule((d,D), B), % pick a rule
	append(B, RestGs, NewGs), % FIXME: solve constraints first?
	NextState = (NewGs, As, Cs, Ns, Q, Info).

h_solve_one(b, gpr(_A,_B), RestGs, As, Cs, Ns, Q, Info, NextState) :-
        !,
        checkAs(As),
%        learning_score_c(([gpr(A,B) | RestGs], As, Cs, Ns, Q, Info), (RestGs, As, Cs, Ns, Q, Info), SCORE),
%        set_info(Info, Info2, scr, SCORE),
        NextState = (RestGs, As, Cs, Ns, Q, Info).

h_solve_one(b, minimize(X), RestGs, As, Cs, Ns, Q, Info, NextState) :- !,
        minimize(X), % backtrackable
        NextState = (RestGs, As, Cs, Ns, Q, Info).

h_solve_one(b, maximize(X), RestGs, As, Cs, Ns, Q, Info, NextState) :- !,
        minimize(X), % backtrackable
        NextState = (RestGs, As, Cs, Ns, Q, Info).

h_solve_one(b, B, RestGs, As, Cs, Ns, Q, Info, NextState) :-
	call(B), % backtrackable
	NextState = (RestGs, As, Cs, Ns, Q, Info).

h_solve_one(e, X = Y, RestGs, As, Cs, Ns, Q, Info, NextState) :-
	call(X = Y),
	NextState = (RestGs, As, Cs, Ns, Q, Info).

h_solve_one(i, X=/=Y, RestGs, As, Cs, Ns, Q, Info, NextState) :-
	propagate_inequalities([X=/=Y], Cs, NewCs), 
	NextState = (RestGs, As, NewCs, Ns, Q, Info).

h_solve_one(n, \+ G, RestGs, As, Cs, Ns, Q, Info, NextState) :-
	h_solve_one(-, fail([], [G]), RestGs, As, Cs, Ns, Q, Info, NextState).

h_solve_one(f, FC, RestGs, As, Cs, Ns, Q, Info, NextState) :-
	propagate_finite_domain_constraints([FC], Cs, NewCs),
	NextState = (RestGs, As, NewCs, Ns, Q, Info).

h_solve_one(r, RC, RestGs, As, Cs, Ns, Q, Info, NextState) :-
	propagate_real_domain_constraints([RC], Cs, NewCs), 
	NextState = (RestGs, As, NewCs, Ns, Q, Info).

h_solve_one(-, fail(Vs, Lits), RestGs, As, Cs, Ns, Q, Info, NextState) :-
	(safe_select_failure_literal(Lits, Vs, (Type, L), RestLits) ->
		h_fail_one(Type, L, Vs, RestLits, RestGs, As, Cs, Ns, Q, Info, NextState)
	;
		(Lits == [] ->
			fail % fail normally
		;
			% if "Lits" contains only finite domain constraints and/or 
			% real constraints, then we can try to test for satisfiability.
			% if it is satisfiable, then it fails normally (i.e. falsity 
			% can be derived); otherwise if it is not satisfiable, then
			% it succeeds.
			(partition_failure_literals(Lits, FCs, RCs) ->
				(\+ \+ (
							% can be satisfied?
							propagate_finite_domain_constraints(FCs, Cs, Cs1),
							propagate_real_domain_constraints(RCs, Cs1, _)
						) ->
					% yes, satisfiable.  so fail
					fail
				;
					% no, and hence no floundering, so continue the reasoning
					NextState = (RestGs, As, Cs, Ns, Q, Info)
				)
			; % indeed, it flounders.
				write('Floundered: '), write(fail(Vs, Lits)), nl,
				fail
			)
		)
	).

h_fail_one(a, arwr(_AB), Vs, RestLits, RestGs, As, Cs, Ns, Q, Info, NextState) :- %log((a, arwr(AB), Vs, RestLits, RestGs, As, Cs, Ns, Q, Info, NextState)), 
        !,
                                                                                 %log(resolve_failure_abducible_with_delta(As, arwr(AB), Vs, RestLits, FailureGoals)),
        %resolve_failure_abducible_with_delta(As, arwr(AB), Vs, RestLits, FailureGoals),
                                %log(append(FailureGoals, RestGs, NewGs)),
        %append(FailureGoals, RestGs, NewGs),
        NewNs = Ns,
        NextState = ([(-, fail(Vs, RestLits)) | RestGs], As, Cs, NewNs, Q, Info).


h_fail_one(a, A, Vs, RestLits, RestGs, As, Cs, Ns, Q, Info, NextState) :-
	resolve_failure_abducible_with_delta(As, A, Vs, RestLits, FailureGoals),
	append(FailureGoals, RestGs, NewGs),
	NewNs = [fail(Vs, [(a, A)|RestLits])|Ns],
	NextState = (NewGs, As, Cs, NewNs, Q, Info).

h_fail_one(d, D, Vs, RestLits, RestGs, As, Cs, Ns, Q, Info, NextState) :-
	findall(H-B, (rule((d, H), B), unifiable(H, D)), Rules),
	resolve_failure_non_abducible_with_rules(Rules, D, Vs, RestLits, FailureGoals),
	append(FailureGoals, RestGs, NewGs),
	NextState = (NewGs, As, Cs, Ns, Q, Info).

h_fail_one(b, B, Vs, RestLits, RestGs, As, Cs, Ns, Q, Info, NextState) :-
	findall(B-[], call(B), Rules),
	resolve_failure_non_abducible_with_rules(Rules, B, Vs, RestLits, FailureGoals),
	append(FailureGoals, RestGs, NewGs),
	NextState = (NewGs, As, Cs, Ns, Q, Info).

h_fail_one(e, X = Y, Vs, RestLits, RestGs, As, Cs, Ns, Q, Info, NextState) :-
	((var(X), strictmember(Vs, X)) ->
		% X is uni. quant. and we don't care about Y
		strictdelete(Vs, X, NewVs),
		call(X = Y),
		h_solve_one(-, fail(NewVs, RestLits), RestGs, As, Cs, Ns, Q, Info, NextState)
	; (var(Y), strictmember(Vs, Y)) ->
		% Y is uni. quant. but X is not
		strictdelete(Vs, Y, NewVs),
		call(Y = X),
		h_solve_one(-, fail(NewVs, RestLits), RestGs, As, Cs, Ns, Q, Info, NextState)
	; var(X) ->
		% X is ex. quant. 
		% NB: by the safe selection strategy, "Y" doesn't contain any
		% universally quantified variable
		(
			% try to succeed in the inequality
			propagate_inequalities([X=/=Y], Cs, NewCs),
			NextState = (RestGs, As, NewCs, Ns, Q, Info)
		;
			% try to succeed in the equality and fail one of the rest literals
			call(X = Y),
			h_solve_one(-, fail(Vs, RestLits), RestGs, As, Cs, Ns, Q, Info, NextState)
		)
	; var(Y) ->
		% Y is ex. quant. but X is not
		% NB: by the safe selection strategy, "X" doesn't contain any
		% universally quantified variable
		(
			% try to succeed in the inequality
			propagate_inequalities([Y=/=X], Cs, NewCs),
			NextState = (RestGs, As, NewCs, Ns, Q, Info)
		;
			% try to succeed in the equality and fail one of the rest literals
			call(Y = X),
			h_solve_one(-, fail(Vs, RestLits), RestGs, As, Cs, Ns, Q, Info, NextState)
		)
	; % one of them is a variable
		(unifiable(X, Y, Es) ->
			maplist(wrap_atom, Es, EsW),
			append(EsW, RestLits, NewBody),
			h_solve_one(-, fail(Vs, NewBody), RestGs, As, Cs, Ns, Q, Info, NextState)
		;
			% this literal fails trivially
			NextState = (RestGs, As, Cs, Ns, Q, Info)
		)
	).

h_fail_one(i, X=/=Y, Vs, RestLits, RestGs, As, Cs, Ns, Q, Info, NextState) :-
	(
		call(X = Y),
		NextState = (RestGs, As, Cs, Ns, Q, Info)
	;
		propagate_inequalities([X=/=Y], Cs, NewCs),
		h_solve_one(-, fail(Vs, RestLits), RestGs, As, NewCs, Ns, Q, Info, NextState)
	).

h_fail_one(n, \+ G, Vs, RestLits, RestGs, As, Cs, Ns, Q, Info, NextState) :-
	(
		% Case 1:
		NextState = ([G|RestGs], As, Cs, Ns, Q, Info)
	;
		% case 2:
                RestLits \= [], %DOMENICO-5-Sept
		G = (Type, L),
		h_fail_one(Type, L, [], [], [(-, fail(Vs, RestLits))|RestGs], As, Cs, Ns, Q, Info, NextState)
	).

h_fail_one(f, C, Vs, RestLits, RestGs, As, Cs, Ns, Q, Info, NextState) :-
	fd_entailment(C, B), % reitification
	(B == 0 -> 
		% the constraint can never hold anyway
		NextState = (RestGs, As, Cs, Ns, Q, Info)
	; B == 1 ->
		% the constraint always hold
		h_solve_one(-, fail(Vs, RestLits), RestGs, As, Cs, Ns, Q, Info, NextState)
	; % the constraint can hold and can unhold
		% then we either (a) fail it or (b) succeed it and fail the rest
		(
			negate_finite_domain_constraint(C, NC), 
			propagate_finite_domain_constraints([NC], Cs, NewCs),
			NextState = (RestGs, As, NewCs, Ns, Q, Info)
		;
			propagate_finite_domain_constraints([C], Cs, NewCs),
			h_solve_one(-, fail(Vs, RestLits), RestGs, As, NewCs, Ns, Q, Info, NextState)
		)
	).

h_fail_one(r, C, Vs, RestLits, RestGs, As, Cs, Ns, Q, Info, NextState) :-
	negate_real_domain_constraint(C, NC),
	C = {Cons}, NC = {NCons},
	(entailed(NCons) -> 
		% the constraint can never hold anyway
		NextState = (RestGs, As, Cs, Ns, Q, Info)
	; entailed(Cons) ->
		% the constraint always hold
		h_solve_one(-, fail(Vs, RestLits), RestGs, As, Cs, Ns, Q, Info, NextState)
	; % we either (a) fail it, or (b) succeed it and fail the rest.
		(
			propagate_real_domain_constraints([NC], Cs, NewCs),
			NextState = (RestGs, As, NewCs, Ns, Q, Info)
		;
			propagate_real_domain_constraints([C], Cs, NewCs),
			h_solve_one(-, fail(Vs, RestLits), RestGs, As, NewCs, Ns, Q, Info, NextState)
		)
	).

% --- Auxiliary Predicates for the Meta-Interpreters ---

% initialise(Query, InitGs, Abducibles, Constraints, Denials, Info)
initialise(Query, InitGs, [], ([], [], []), InitNs, Info2) :-
	maplist(wrap_literal, Query, QueryW),
	collect_static_ics(FGs, InitNs),
	append(QueryW, FGs, InitGs),
	debug_hook(init_state(Query), Info),
        set_info(Info, Info2, depth, 0).

collect_static_ics(FailureGoals, DynamicDenials) :-
	findall(IC, ic(IC), ICs),
	partition_ics(ICs, FailureGoals, DynamicDenials).

partition_ics([], [], []).
partition_ics([IC|Rest], AllFs, AllNs) :-
	partition_ics(Rest, Fs, Ns),
	term_variables(IC, Vs),
	(select((a, A), IC, RestBody) ->
		% contains an abducible, so move it to the front to form a dynamic denial
		AllNs = [fail(Vs, [(a,A)|RestBody])|Ns],
		AllFs = Fs
	;
		% otherwise, keep it as a failure goal
		AllFs = [(-, fail(Vs, IC))|Fs],
		AllNs = Ns
	).

resolve_abducible_with_delta([], _, []).
resolve_abducible_with_delta([H|T], A, [(A=/=H)|IEs]) :-
	unifiable(H, A), !,
	resolve_abducible_with_delta(T, A, IEs).
resolve_abducible_with_delta([_|T], A, IEs) :-
	resolve_abducible_with_delta(T, A, IEs).

resolve_abducible_with_dynamic_ics([], _, []).
resolve_abducible_with_dynamic_ics([fail(Vs, [(a, Left)|Rest])|RestNs], A, [F|Fs]) :-
	unifiable(Left, A), !,
	rename_universal_variables(Vs, Left-Rest, RnVs, RnLeft-RnRest),
	unifiable(RnLeft, A, Bindings),
	bind_universal_variables(Bindings, RnVs, RemainingBindings, RemainingRnVs),
	maplist(wrap_atom, RemainingBindings, Es),
	append(Es, RnRest, NewRest),
	F = (-, fail(RemainingRnVs, NewRest)),
	resolve_abducible_with_dynamic_ics(RestNs, A, Fs).
resolve_abducible_with_dynamic_ics([_|RestNs], A, Fs) :-
	resolve_abducible_with_dynamic_ics(RestNs, A, Fs).

bind_universal_variables([], Vs, [], Vs).
bind_universal_variables([(A=B)|T], Vs, NewEs, NewVs) :-
	(strictmember(Vs, A) ->
		% A is universally quantified, so bind it
		strictdelete(Vs, A, Vs1),
		call(A=B),
		bind_universal_variables(T, Vs1, NewEs, NewVs)
	;
		% A is not universally quantified, so keep the equality
		bind_universal_variables(T, Vs, Es1, NewVs),
		NewEs = [(A=B)|Es1]
	).

% rename all the universally quantified variables in a given term.
% HACK: the term_variables/2 in YAP and that in SICStus behave
% differently.  Therefore, the implementation of "rename_universal_variables" 
% can't rely on them any more.  Below is a new implementation that should work
% on both systems

:- dynamic mycopy/1.
rename_universal_variables(VsIn, FormIn, VsOut, FormOut) :-
	assert(mycopy(VsIn-FormIn)), retract(mycopy(VsOut-FormOut)),
	rebind_existential_variables(FormIn, VsIn, FormOut).
rebind_existential_variables(Pin, VsIn, Pout) :-
	(atomic(Pin) ->
		true % nothing to bind
	; var(Pin) ->
		(strictmember(VsIn, Pin) ->
			% this is a universally quantified variable. so do not bind it
			true
		;
			% this is an existentially quantified variable, so BIND it.
			Pout = Pin
		)
	; % this is a compound term
		Pin =.. [F|ArgsIn],
		Pout =.. [F|ArgsOut],
		rebind_existential_variables_recursively(ArgsIn, VsIn, ArgsOut)
	).
rebind_existential_variables_recursively([], _, []).
rebind_existential_variables_recursively([H1|T1], VsIn, [H2|T2]) :-
	rebind_existential_variables(H1, VsIn, H2),
	rebind_existential_variables_recursively(T1, VsIn, T2).

safe_select_failure_literal(Lits, Vs, (Type, L), RestLits) :-
	select((Type, L), Lits, RestLits), 
	safe_failure_literal(Type, L, Vs), !.

% 17/09/2010 Bug fixed: thanks Domenico Corapi
partition_failure_literals([], [], []).
partition_failure_literals([(f,C)|T], [C|FCs], RCs) :-
	partition_failure_literals(T, FCs, RCs).
partition_failure_literals([(r,C)|T], FCs, [C|RCs]) :-
	partition_failure_literals(T, FCs, RCs).

safe_failure_literal(n, L, Vs) :-
	% negative literal must not contain universally quantified 
	% variables
	variables_within_term(Vs, L, []).
safe_failure_literal(i, L, Vs) :-
	% simliary inequality must not involve universally quantified
	% variables
	variables_within_term(Vs, L, []).
safe_failure_literal(f, L, Vs) :-
	% finite domain constraint with universally quantified variables
	% should not be selected at first
	variables_within_term(Vs, L, []).
safe_failure_literal(r, L, Vs) :-
	% similarly real domain constraint with universally quantified 
	% variables should not be selected at first
	variables_within_term(Vs, L, []).
safe_failure_literal(e, X = Y, Vs) :-
	% for equalitiy it is a bit trickier
	% FIXME (low priority): at the moment, the equality store does
	% not handle inequalities with universally quantified variables.
	% Thus, equalities between one existentially quantified variable
	% and a compound term with universally quantified variables should
	% not be selected.
	\+ (
		% the following means "X is an existentially quantified variable
		% but Y is a compound term with universally quantified variables"
		var(X), \+ strictmember(Vs, X),
		\+ var(Y), \+ variables_within_term(Vs, Y, [])
	), 
	\+ (
		% the following means "Y is an existentially quantified variable
		% but X is a compound term with universally quantified variables"
		var(Y), \+ strictmember(Vs, Y),
		\+ var(X), \+ variables_within_term(Vs, X, [])
	).
% finally, everything else is safe.
safe_failure_literal(a, _, _).
safe_failure_literal(b, _, _).
safe_failure_literal(d, _, _).

resolve_failure_abducible_with_delta([], _, _, _, []).
resolve_failure_abducible_with_delta([H|T], A, Vs, RestLits, [FailureGoal|FailureGoals]) :-
	unifiable(H, A), !,
	rename_universal_variables(Vs, A-RestLits, RnVs, RnA-RnRestLits),
	unifiable(RnA, H, Es),
	bind_universal_variables(Es, RnVs, RemainingEs, RemainingRnVs),
	maplist(wrap_atom, RemainingEs, EsW),
	append(EsW, RnRestLits, NewBody),
	FailureGoal = (-, fail(RemainingRnVs, NewBody)),
	resolve_failure_abducible_with_delta(T, A, Vs, RestLits, FailureGoals).
resolve_failure_abducible_with_delta([_|T], A, Vs, RestLits, FailureGoals) :-
	resolve_failure_abducible_with_delta(T, A, Vs, RestLits, FailureGoals).

resolve_failure_non_abducible_with_rules([], _, _, _, []).
resolve_failure_non_abducible_with_rules([H-B|T], P, Vs, RestLits, 
		[FailureGoal|FailureGoals]) :-
	rename_universal_variables(Vs, P-RestLits, RnVs, RnP-RnRestLits),
	term_variables(H-B, NewVs),
	append(NewVs, RnVs, AllVs), % no need to use union/3 here as "NewVs" must be fressh
	unifiable(H, RnP, Es),
	bind_universal_variables(Es, AllVs, RemainingEs, RemainingVs),
	maplist(wrap_atom, RemainingEs, EsW),
	append(EsW, B, NewB), append(NewB, RnRestLits, AllBody),
	FailureGoal = (-, fail(RemainingVs, AllBody)),
	resolve_failure_non_abducible_with_rules(T, P, Vs, RestLits, FailureGoals).

% --- Utilities for Inequality and Constraint stores ---

propagate_inequalities(Es, (E, F, R), (NewE, NewF, NewR)) :-
	add_inequalities(Es, (E, F, R), (NewE, NewF, NewR)).

add_inequalities([], Cs, Cs).
add_inequalities([(X=/=Y)|T], (Ei, Fi, Ri), (Eo, Fo, Ro)) :-
	(finite_domain_inequality(X=/=Y) ->
		% push it to the finite domain constraint store
		fd_entailment((X #\= Y), B),
		(B == 1 ->
			% already reitified, so can discard it
			add_inequalities(T, (Ei, Fi, Ri), (Eo, Fo, Ro))
		; % otherwise, try to add it to the finite domain constraint store
			call(X #\= Y),
			add_inequalities(T, (Ei, [(X #\= Y)|Fi], Ri), (Eo, Fo, Ro))
		)
	; real_domain_inequality(X=/=Y) ->
		% push it to the real domain constraint store
		(entailed(X =\= Y) ->
			% already entailed, so can discard it
			add_inequalities(T, (Ei, Fi, Ri), (Eo, Fo, Ro))
		; % otherwise, try to add it to the real domain constraint store
			call({X =\= Y}),
			add_inequalities(T, (Ei, Fi, [{X =\= Y}|Ri]), (Eo, Fo, Ro))
		)
	; % a real inequality?
		(unifiable(X, Y) ->
			% still need to keep it as the inequality could be falsified
			call(X=/=Y),
			add_inequalities(T, ([(X=/=Y)|Ei], Fi, Ri), (Eo, Fo, Ro))
		;
			add_inequalities(T, (Ei, Fi, Ri), (Eo, Fo, Ro))
		)
	).

propagate_finite_domain_constraints(Fs, (E,F,R), (NewE,NewF,R)) :-
	add_finite_domain_constraints(Fs, F, F1),
	finite_domain_soundness_check(F1),
	selectlist(nonground, F1, NewF),
	selectlist(nonground, E, NewE).

add_finite_domain_constraints([], F, F).
add_finite_domain_constraints([H|T], Fin, Fout) :-
	fd_entailment(H, B),
	(B == 1 ->
		% can discard it
		add_finite_domain_constraints(T, Fin, Fout)
	;
		call(H),
		add_finite_domain_constraints(T, [H|Fin], Fout)
	).

domain_bounded(X) :- fd_size(X, S), S \== sup. 

finite_domain_soundness_check(F1) :-
	term_variables(F1, Vs),
	selectlist(domain_bounded, Vs, BVs),
	(get_value(lbl, true) ->
		labeling([], BVs)
	;
		\+ \+ labeling([], BVs)
	).

propagate_real_domain_constraints(Rs, (E,F,R), (E,F,NewR)) :-
	add_real_domain_constraints(Rs, R, NewR).

add_real_domain_constraints([], R, R).
add_real_domain_constraints([H|T], Rin, Rout) :-
	H = {C},
	(entailed(C) ->
		% can discard it
		add_real_domain_constraints(T, Rin, Rout)
	;
		call(H),
		add_real_domain_constraints(T, [H|Rin], Rout)
	).

negate_finite_domain_constraint(C1 #\/ C2, NC1 #/\ NC2) :-
	negate_finite_domain_constraint(C1, NC1),
	negate_finite_domain_constraint(C2, NC2).
negate_finite_domain_constraint(C1 #/\ C2, NC1 #\/ NC2) :-
	negate_finite_domain_constraint(C1, NC1),
	negate_finite_domain_constraint(C2, NC2).
negate_finite_domain_constraint(#\ C, C).
negate_finite_domain_constraint(X #< Y, X #>= Y).
negate_finite_domain_constraint(X #> Y, X #=< Y).
negate_finite_domain_constraint(X #=< Y, X #> Y).
negate_finite_domain_constraint(X #>= Y, X #< Y).
negate_finite_domain_constraint(X #\= Y, X #= Y).

negate_real_domain_constraint({C1 , C2}, {NC1 ; NC2}) :-
	negate_real_domain_constraint({C1}, {NC1}),
	negate_real_domain_constraint({C2}, {NC2}).
negate_real_domain_constraint({C1 ; C2}, {NC1 , NC2}) :-
	negate_real_domain_constraint({C1}, {NC1}),
	negate_real_domain_constraint({C2}, {NC2}).
negate_real_domain_constraint({X < Y}, {X >= Y}).
negate_real_domain_constraint({X > Y}, {X =< Y}).
negate_real_domain_constraint({X =< Y}, {X > Y}).
negate_real_domain_constraint({X >= Y}, {X < Y}).
negate_real_domain_constraint({X =\= Y}, {X = Y}).
negate_real_domain_constraint({X = Y}, {X =\= Y}).
negate_real_domain_constraint({X =:= Y}, {X =\= Y}).

extract_constraints((Es, Fs, Rs), Cs) :-
	% collect inequalities
	term_variables(Es, EsVs),
	maplist(inequalities, EsVs, Ess),
	flatten_lists(Ess, Es1),
	list_to_ord_set(Es1, Es2),
	% collect ground finite domain constraints
	selectlist(nonground, Fs, Fs1),
	% collect real domain constraints
	term_variables(Rs, RsVs),
	dump(RsVs, RsVs, Rs1),
	append(Fs1, Rs1, Cs1),
	append(Es2, Cs1, Cs).

flatten_lists([], []).
flatten_lists([H|T], L) :-
        flatten_lists(T, L1),
        append(H, L1, L).

